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Abstract 

We report several important observations that underscore the distinctions 
between the constrained-path Monte Carlo method and the continuum and 
lattice versions of the fixed-node method. The main distinctions stem from 
the differences in the state space in which the random walk occurs and in 
the manner in which the random walkers are constrained. One consequence 
is that in the constrained-path method the so-called mixed estimator for the 
energy is not an upper bound to the exact energy as previously claimed. 
Several ways of producing an energy upper bound are given, and relevant 
methodological aspects are illustrated with simple examples. 
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I. INTRODUCTION 



It is arguable that the fixed-nodeS and constrained-pathi quantum Monte Carlo meth- 



ods are the two most powerful and useful simulation techniques for computing accurate 
ground-state (T = K) properties of large systems of interacting quantum particles. As the 
significantly older method, the fixed-node method has been well studied, and its properties 
are well documented. i Less is known about the constrained-path method, but because of our 
recent use of this methodJHl we now can report several important experiences and obser- 
vations that underscore features distinguishing it from both the continuumi!0 and latticed 
versions of the fixed-node method. 

There is a very strong analogy between the fixed-node and constrained-path methods. 
Both, in a sense, are auxiliary-field methods, both project the ground-state wave function 
from a trial wave function by an importance-sampled, branched random walk, and both place 
a constraint on this random walk to prevent the fermion sign problem from rapidly producing 
exponentially growing variances. A number of technical details for their implementation are 
the same. In fact, the formal development of the constrained-path method! relied on the 
existence of the fixed-node method. The three principal differences between the methods 
are (a) the state space where random walks have their support, (b) the manner by which 
random walkers are constrained, and (c) the part of the imaginary-time propagator that is 
stochastically sampled. 

The continuum version of the fixed- node method works in a first quantized representation 
and operates in coordinate space. The basis states are the complete orthonormal set of 
the particle configurations. The constrained-path method works in a second quantized 
representation and operates in Fock space. Its basis states are the over-complete non- 
orthonormal set of Slater determinants 




(1) 



where the a\ = J2j=i &ij c ] creates a fermion in a quasi particle state i defined relative to M 
possible single particle states j created by the operator c], and |0) represents the vacuum. 
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(In the one-band Hubbard model M will be the number of lattice sites.) In this basis 
classes of many-electron wave functions like the BCS wave function are more easily used, 
and many-particle expectation values like superconducting pairing correlation functions are 
more easily evaluated than possible with the fixed-node method.i The ease in the evaluation 
of ground-state observables, for example, is a consequence of the ease in evaluating single- 
particle propagators and using Wick's theoremEi to express any multi-particle propagator 
as a linear combination of products of one-particle propagators. 

The differences in bases generate a difference in the way the random walks are con- 
strained. Both methods rely on a trial state \^t) t° perform the constraint. In the contin- 
uum version of the fixed-node method with r, representing a particle's position, the random 
walks are confined within a surface defined by (R|\I/y) = \I/t(R) = \I/T( r i, r 2, • • • , r 7v) > 0, 
whereas in the constrained-path method, only random walkers \4>) satisfying > are 

permitted. The fixed-node method solves Schrodinger's equation for the ground-state wave 
function inside the nodal surface. Unless that surface is exact, only an approximate solution 
is obtainable. The constrained-path condition, as we will discuss, has different implications. 
In certain cases, including some simple examples detailed below, the constraint is never 
invoked and hence the constrained-path method can sometimes produce the exact solution 
even for systems of interacting fermions. We also give examples where the solution, though 
approximate, is extremely accurate. These examples include a closed-shell Hubbard model 
with a large positive U. 

The resulting stochastic dynamics in the basis space (both coordinate and Slater deter- 
minantal manifolds) is a Markov process generated by a conditional probability connected 
to the imaginary-time r propagator exp[— rif], where H = T + V is the Hamiltonian repre- 
senting the system and, as usual, T and V are the kinetic and potential energy operators, 
respectively. The kinetic energy propagator is non-diagonal in the coordinate basis repre- 
sentation and its action can be viewed as a diffusion process in the basis space. On the other 
hand, in the Slater determinant representation it is the potential energy kernel which, after 
a Hubbard- Stratonovich transformation, generates the Markov chain. As a consequence, 
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non-interacting problems do not suffer from the fermion sign instability. 

The main purpose of the paper is to illustrate that the constrained-path method is not 
equivalent to the fixed-node method in the space of Slater determinants. In Section [TI| we will 



the mixed estimator for the energy and the extent to which it is an upper bound on the 
ground-state energy. In particular, we will conclude and illustrate that in the constrained- 
path method it is not, in general, an upper bound, contrary to previous claims and to the 
fixed-node method. In Section [TV] we provide a correction to the mixed estimator that makes 
it a rigorous upper bound, plus several alternative ways to produce energy estimates that 
are upper bounds. In Section [V|, we conclude by commenting on areas needing additional 
clarification and several other differences between the methods. Some of these differences 
will be illustrated in the Appendix where we present a constrained-path simulation on a 
small toy problem for which many of the details can be generated analytically. Of particular 
emphasis here will be the effects of matrix stabilization. 



Both the fixed- node and constrained-path methods project the ground state |^ ) from 
the long-time solution of the imaginary-time r representation of Schrodinger's equation 
specified by a Hamiltonian H 



summarize the essential mathematical structure of both methods. In Section [Til], we discuss 



II. SUMMARY OF THE TWO METHODS 




(2) 



Provided Nq = (\E'o|^ / (0)) 7^ and H is time-independent, the formal solution 



tf(r)) =e- T{H - Eo) \V(0)) 



(3) 



has the property 




(4) 
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On the computer this limit is accomplished iteratively 

where At is a small number, nAr = r, n is the current number of iterations (often called 
time slices), and Et is a trial guess at the ground-state energy E Q . The iterative process 
is converted into a stochastic sampling process. As the matrix elements of the propagator 
exp[— t(H — Et)] between different anti-symmetric wave functions are not always positive 
definite, constraints in the sampling are necessary to insure this. Additionally, importance 
sampling is also required to control the variances of computed results. If Et is adjusted so 
that it equals E , then as r — > oo, the iteration becomes stationary, i.e. d\^>)/dr = and 

For simplicity, we will exclude branching from our discussions and consider only time- 
reversal symmetric Hamiltonians, that is, real symmetric operators for which the ground- 
state wave functions can always be chosen to be real. This analysis leaves out the very 
important case of systems in the presence of external magnetic fields.ll! Here we compare the 
constrained-path method to the continuum fixed-node approach. There are some technical 
differences between the continuum fixed-node methodi'0 and the lattice version0 we prefer 
to omit. In this regard we comment that the constrained-path method does not distinguish 
between lattice and continuum fermions: both are treated on an equal footing. 



A. Fixed-Node method 

In the fixed-node method, one represents the ground state as |^ ) = J2n |R)(R|^o) = 
J2n ^o(R-) |R) where \l/o(R) > 0. Asymptotically, the Monte Carlo procedure samples from 
the distribution P(R) = # (R)/ Er *o(R)- 

In the fixed-node method, one projects the iteration onto the basis of particle configura- 
tions {|R)} 



5 



LA-UR-98-5289 February 1, 2008 submitted to Phys. Rev. B 

Projecting this equation onto (R'| and inserting J^n |R) (R| = 1 leads to 

(R|^'> = #'(R') = ]T (R'|e- Ar{H - BT) |R)^(R) , (7) 

R 

and correspondingly the imaginary-time Schrodinger's equation becomes 

" d ^f T ,T) = t-^V 2 + V(R) - Et]*(R, t) (8) 

where D = h 2 /2m, m is the fermion mass, and V^(R) is the potential energy. 

^(R) must be positive to be interpreted as the limiting probability distribution of the 
Markov chain. Fixing the node forces this by prohibiting any change in the particle config- 
uration R — > R' that changes the sign of ^(R). We will denote the ground state produced 
under the constraint, i.e. under the fixed-node condition, as |\l/ c ) and the eigenvalue of the 
propagation (Eq. |6[) as exp[— Ar(E g — Et)\- This eigenvalue defines E g , the growth energy. 

After the so-called short-time approximation is made on the kernel of the integral, which 
is equivalent to making a Trotter approximation and a Hubbard-Stratonovich transformation 
on the exponential of the kinetic energy! the positivity of (R'| exp[— Ar(H — Et)]\R) is 
trivially satisfied so it can be interpreted as a transition probability defining a Markov 
chain. We will call the resulting approximation K(R — > R'). 

It is critical to importance sample in order to reduce statistical fluctuations, especially 
when the potential l^(R) has some singularities, for example, like the 1/r Coulomb singu- 
larity. This means we generate a new distribution \I/ C (R) = \l/r(R)\l/ c (R) satisfying 

*' C (R') = £ K(R -> R')* C (R) , (9) 

R 

where K(R — > R') = ^(R^i^R — > R')/\I/t(R). The new configurations are now sampled 
with a different probability. The new distribution also satisfies a different equation of motion 

<# c (R,r) - 



dr 



DW C (R, t) + DV- [* c (R, r)F(R)] + (E L (R) - E T )V C (R, r) 



(10) 

where the "quantum drift" F = 2Vln^ T and the "local energy" E L (R) = ff* r (R)/* r (R). 
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The Monte Carlo procedure represents the multi-dimensional integral as a set of random 
walkers {|R)} where each member of the set is a different allowed particle configuration. A 
new configuration |R') is sampled from K(R, — > R') and rejected, thereby terminating this 
random walker, if the resulting value of the wave function is negative. 

Since f c (R, r) = \& C (R, r)/\E' r (R), a variational upper bound to the true energy is 



At large times walkers are distributed with a probability density ^t{R)^c{R), and both ty? 
and \l/ c go to zero linearly near the nodal surface. Since the Hamiltonian and the constraint 
are all local operators, this implies that the growth energy E g is equal to the mixed estimate 
of the energy E n ^ 1 " 



The constrained propagator is identical to the exact one except near the nodal surface where 
the constraint acts. The constraint discards contributions that are orthogonal to both 
and \l/ c , and hence this region gives no net contribution to either {^t\H\^ c ) or (^ c \H\ty c ). 
Therefore, the variational estimate of the energy E v is identical to E g and E m , and all are 
variational upper bounds. E m is more easily, accurately, and efficiently computed than E g 
or E v . 

Several characteristics of the fixed- node method are: 1) The nodal surface of \I/ C (R) is 
exactly the same as that of ^y(R); 2) The exact ground-state energy is obtained only if the 
nodal surface of \I/t(R) is exact; and 3) even for the trivial case of V^(R) = 0, unless the 
exact nodal surface is used, only an approximate solution is produced. 

B. Constrained-Path Method 

In the constrained-path method, one represents the ground state as |^o) = S^c^l^) 
where the Slater determinants \<f>) are chosen so that all > 0. Asymptotically, the Monte 
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Carlo procedure samples from the distribution 7r(0) = c^j J2</> c <t>- The decomposition of |\I/o) 
in terms of the |0}'s is not unique. One could just a well have |^ ) = 12$ d^cj)) where > 0. 
We will simply write |^o) = 12$ !</>)■ 

The constrained-path method works in a basis of Slater determinants. Again one iterates 
Eq. |6], placing constraints on the random walks. A different kind of constraint is needed 
because a different basis is used. Here the sign problem is caused by transitions from a region 
where the overlap (^t\4>) is positive to a region where it is negative. These two regions are 
not physically distinguishable; they involve merely the exchange of fermions. Hence an 
arbitrary wave function can always be expanded in the restricted bases where (^t\4>) is 
purely positive or purely negative. The original propagation mixes these two degenerate 
bases indiscriminately, causing a sign problem. 

To break this plus-minus symmetry, the random walks are constrained to the region 
(\I/t|0) > 0. This is an approximation because in general a wave function will have both 
positive and negative coefficients c$ when expressed in this basis. However, the constrained 
propagation yields all the c$ > 0. To compare with the fixed-node method, we sketch 
some additional details: After the application of a Trotter approximation and Hubbard- 
Stratonovich transformation, the iterative equation becomes 

|*') = £P(x)5(x)|*), (13) 

X 

where x (the Hubbard- Stratonovich field) is to be interpreted as a multi-dimensional random 
variable distributed according to -P(x), and -B(x) is an operator approximating exp[— AtH] 
for a given value of the random variable, whose general structure is a product of exponentials 
of one-body operators. -B(x) has the property of transforming one Slater determinant into 
another. The Monte Carlo method is used to evaluate the multi-dimensional integration by 
using multiple random walkers \<p), and for each walker, sampling a x from -P(x) and then 
generating a new walker 

|0')=£(x)|0> . (14) 
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Thus, if ($ T \4>') > 0, (* T |5(x)|0) > 0. 

Since the basis of Slater determinants is non-orthogonal and over-complete, each member 
of the basis, in general, is a linear combination of the others, i.e., \<p) = Yl's 1 a <t>'\ < f ) ')- A prime 
is on the summation symbol because while the summation may be over an infinite number 
of Slater determinants, the ones used need not exhaust the basis. While \^t) ma Y constrain 
a \4>) to be in the "positive" set, this |0) can overlap with a state in the "negative" set. In 
contrast to the fixed-node condition, the constrained-path condition does not separate the 
basis into orthogonal sets. Whereas the fixed-node condition must produce an approximate 
solution unless the nodes are exact, the constrained-path method can sometimes produce 
the exact solution even if the constraining wave function is approximate and has the wrong 
nodal surface in configuration space. 

Importance sampling is also implemented in the constrained-path method. With |0) = 
(\1/t|0)|0) the iteration on each walker becomes 

\p) = B{x)\$) , (15) 

but now the random variable x is sampled from -P(x) oc (\E't|0 / )-P( x )/(^ / tI0)- We have 

Again a variational estimate of the energy E v can be constructed from \^ c ) 

E -=-4<kr =E '- E °- (16) 

but now the connection among E v , E g and E m is unclear because the constraint discards 
configurations which are orthogonal to \^t) and these discarded configurations are not 
necessarily orthogonal to |\P C ). As we will argue in the next section, the mixed estimator 
is not always an upper bound to E . This retracts previous claims of E m being an upper 
bound.S 

Several characteristics of the constrained-path method are: 1) The nodal surface of (<p\^f c ) 
is not the same as that of (0|\I/t); 2) In some cases, the exact ground-state energy can be 
obtained even if the nodal surface of is approximate; and 3) For the trivial case of 

V(R) = 0, the exact solution is produced. 
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Perhaps the best known examples demonstrating the second point are the half-filled 
positive-f/ Hubbard models and negative-?/ Hubbard models, two classes of models that do 
not have a sign problem. To illustrate the first two characteristics of the CPMC method, 
we consider the following half-filled positive-C/ Hubbard model 



(17) 

<r=U i=l 

which is also a simple model for a Heitler-London molecule. The two-particle ground state 
is given by 



1 ~ / + | j j. \ U — Eq / f f + f \ 

-j====== |t [c n c n + c 2T c 2i J + ^_ [c^c 2[ + c 2T c u J 



|0>, 



(18) 



where t = V2t, and the ground-state energy is E = U/2 — y (U/2) 2 + 2t 2 . This state cannot 
be represented by a single Slater determinant, unless U — 0. (See Eq. 1.) 

Since we want to study the nodal structure of different states, we need to parameterize 
the different iable manifold of Slater determinants of two particles. We choose coordinates 
such that a generic point in the manifold (0 1 ,9 2 ) corresponds to the normalized Slater 
determinant 



|0) = (cos#i cJ T + sin 6^ c\^j (cos6 2 cj^ + sin# 2 4|) |0) • 



(19) 



Alternatively, from Eq. 1, we can represent this state by the product of two 2x1 matrices, 
i.e., <3> = <3> T <3>! where 



<3>t = 



and 



Then, 



l*o> 



fit 2 + {U- E y 



( a ^ 

COS V\ 

sin 9i 



^ cos 62 ^ 
sin 9o 



t cos(^i - 9 2 ) + 



(20) 



(21) 



' E ° sin(^ + 9 2 ) 



V2 



(22) 
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In Fig. |l| we display contour plots of this function for different values of U. Clearly the nodal 
surfaces of (<p\^o) are different for the various values of U. Nevertheless, in the absence of 
importance sampling, one can prove analytically that for any U, (0|^t>) = (<f)\^o(U = 0)) 
remains positive during the whole imaginary-time evolution; that is, the nodal constraint is 
never invoked, and therefore, the exact solution is obtained after a large-r projection with 
the result that the nodal surfaces of (<p\^ c ) and {4>\^o) are the same but different from 
(0|*t)- 

III. MIXED ESTIMATOR OF THE ENERGY 

Independent of a quantum Monte Carlo process, when is the mixed estimator for the 
energy (Eq. [12]) an upper bound for the ground-state energy or even the exact value? Three 
cases are apparent: 

1. If |* r ) = |* c ), then E m > E . 

2. If |* r ) ee |tf ), or |* c ) ee \V ), then E m = E . 

3. If \^f c ) = U 2ti |\I/t); where n is an integer, U is a Hermitian operator, and [H, U] = 0, 
then E m > E . 

Case 1 is simply the Rayleigh-Ritz variational principle. Case 2 is perhaps the most impor- 
tant feature of the mixed estimator: a good approximation to the ground-state wave function 
will produce a good approximation to the ground-state energy. Case 3 is what happens in 
quantum Monte Carlo simulations: in principle, with U = exp[—ArH], an upper bound 
on the energy is automatically produced. In the simulations, U 2n — > U 2n • • • U2U1 where 
Uj = exp[— At Hi] with Hi representing an effective Hamiltonian satisfying [H, U»] 7^ 0, so 
in general E m is not a rigorous bound. But since [H, Uj] ~ 0, E m is in general expected to 
be a good estimate and a bound. Clearly to the extent that the constrained-path method 
in principle can produce the exact state vector, the mixed estimate of energy can be exact. 
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On general grounds we can say that if our constrained evolution defines a Markov process 
with a stationary distribution |^ c ), such that 



and H eS = H + 5H, then 



H eS \* c ) = E g \* c ) , (23) 



' =E m + 5E m >E , (24) 



where 

^ m " • (25) 

It is clear that if 5E m < 0, E rn is an upper bound to the ground-state energy. However, in 
general, this is not necessarily the case. 

It is interesting to mention that in the usual fixed-node approach, where the state space 
manifold is the coordinate space, SH represents a hard-wall potential, i.e. it is infinite on 
the set of configurations {|Ry)} defined by (Rt\^t) = 0. Then E m is an upper bound to 
Eq. We can see this by minimizing the following constrained functional 

\^ c ),(^c\;vB. T ,vk 7 

(9 e \H\* e ) - E(V C \V C ) - E^R t ( R tI*c> - ]>>r t <* c |Rt> , (26) 

where (Rt\^t) = ^(Rt) = defines the nodal surface Mr- The resulting Euler's equations 
are 

#|* c ) = £|* c )+5> Rt |R t ) , (27) 



R 



T 



(y c \H=(y c \E + Y,vk T ( R T\ , (28) 



R 



T 



which lead to 



* c (Rr) = , (29) 



'-WW-'" (30) 
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From the first equation 



(R|#|* c > = E (R|* c )+$> Rr (R|R T ) , 



(31) 



£(R|#|R')(R|* C ) = E * C (R) + 5> Rt 5 rr , 



■T 1 



(32) 




H(R) * C (R) = E * C (R) + W] 



RR t • 



(33) 



Solving the constrained (fixed-node) problem is equivalent to solving H(FL) \I/ C (R) = 
E \I/ C (R) within the region where \I/*r(R) has a definite sign, with the boundary condi- 
tion ^/ c (Rj') = 0. In this way \I/ C (R) is a continuous function of R with discontinuous 
derivative at R = Ry. 

If we try to minimize a similar functional, but we use the representation of Slater deter- 
minants then an extremely non-local term, which is not easily handled, appears in the 
resulting Euler equation 



with \I/t[0t] = (0t|^t) = but in general (4>\4>t) ^ S<j><f> T - As before, one can easily prove 
that E m > Eq. In other words, if we had used the exact equivalent of the fixed-node 
constraint, we would have gotten a variational upper bound using E m . It is important to 
stress that the constrained-path condition is a kind of global constraint as opposed to the 
local one that represents the fixed-node constraint in that the constrained-path condition 
does not impose on ((f)\^f c ) the same nodal hypersurface as (4>\^t)- In fact, we have numerical 
examples where (4>\^>t) does not define the exact nodal structure, nevertheless we get the 
exact ground-state energy for H, i.e. \^ c ) = \^o). (See, for instance, the example shown at 
the end of Section 0.) 

IV. ENERGY ESTIMATORS BOUNDING THE GROUND-STATE ENERGY 



£(0|#|0'> * c [0'] = E * c [0'] +Y,^tWt) , 



(34) 
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A. Energy Bounds 

It is possible to construct a variety of other estimators that produce upper bounds to the 
ground-state energy. In the following we assume we have a state |\I/ C ) that is an eigenstate 
of the constrained propagator: 



tf c ) = lim e~ TH |* T ) (35) 



with eigenvalue exp[— rE g 



e -ArH |^ = e -ArE g |^ _ ^ 



The arrow indicates the direction of propagation with the constraint applied to the wave 
function. Only those auxiliary fields that retain a positive overlap with the trial function 
are retained in the sampling. The CPMC paths are not reversible in the standard sense, 
and hence the "time arrow" of the path is significant. In contrast, we denote the original 
unconstrained propagator as exp[— rif], and since H is Hermitian, the full propagation is 
reversible, at least when averaged over paths. The effect of the constraint is simply given 
by the difference between exp[— tH] and exp[— tH]. 
The standard variational upper bound is given by 

^-i^r- Eo ' (37) 

where the function (\P C | is the dual state of |\I/ C ); that is, the constraint is applied in the 
opposite r direction. It is possible to calculate E v directly, for example, by propagating two 
populations of random walkers. These two populations can be used as independent samples 
of ($f c \ and of |\& c ). Since these walkers should be independently evaluated (at least prior to 
the introduction of importance sampling), we label them as ($i c \ and \^ rc )- The importance 
function will presumably have to be a function of all the relevant overlaps, 

/ = I(\(% c \* rc )\, (tt r |*rc), (*ic|*r)) • (38) 
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The overlap of the left and right wave functions may be negative, so we have to assume the 
importance function is only a function of the magnitude of that overlap. In the absence of 
importance sampling, the denominator in Eq. is the sum of the overlap between these 
two wave functions. Hence this term should be large in the importance function. 
One can also evaluate the energy difference E d = E v - E g , given by: 



e -ArE g _ e -ArE v ^ Ar ^ 



(39) 



(*c|*c) 

The numerator in this expression is the result of the constraint: it is simply the overlap 
of (\l/ c | with the state \^d) representing the difference between the full and constrained 



propagation: 



e -ArH _ e -ArH 



l*c> • (40) 



This difference is simply the set of the configurations discarded via the constraint. Ed is 
zero if these discarded configurations are, on average, orthogonal to |^ c ). In the fixed- node 
method the configurations thrown away are by definition orthogonal both to \^t) and |^ c ). 
Here, though, our configurations are in general orthogonal only to \^t), and hence the 
variational and mixed estimates of the energy need not be equal. 

It is still true, however, that E m and E g are equal in the limit of zero time step. The 
density of configurations near the surface {^t\4>) goes to zero rapidly so that the surface 
contributions to the constrained propagator do not give a finite contribution to the growth 
estimate of the energy. 

One can evaluate Ed directly. To evaluate the energy difference we again need inde- 
pendent right- and left-hand wave functions \^f rc ) and (^ c |. Dividing the population into 
two independent halves representing the left- and right-hand states, we can evaluate the 
numerator by taking the overlaps of what is discarded (the difference between the full and 
constrained propagators acting on |\& c )) with the independent solution {^ c \- The denomi- 
nator is just the overlap (\l/ c |\l/ c ) of the two solutions. For larger systems it will likely have 
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high statistical errors, but the numerator may be small enough that this does not matter. 
An explicit numerical example is presented in the next section. 

There are several other ways to produce an energy difference Ed = 0. One possibility 
is to introduce a parameter in the trial state and vary it until Ed = 0. Another possibility 
is changing the constraint. For example, we could discard configurations \<p) for which the 
normalized overlap with the trial wave function is less than or equal to some constant a: 



Varying a until the average overlap of the discarded configurations and the constrained 
solution \^f c ) is greater than or equal to zero produces a variational upper bound for the 
energy, since then E g > E v . In this case it is not necessary to evaluate the denominator of 
Eq. |39|. Also, this procedure is exact for an exact constraining state, since in that case we 
could set a = 0. 

We note that for a ^ the mixed estimate E m is not in general equal to the growth 
estimate E g , as there is a finite surface term that contributes to the difference. In fact, the 
difference E g — E m provides a measure of the error introduced by the constraint. Numerical 
examples are provided in the following section. 

This method is general in that it produces a variational upper bound to the energy for 
any Hamiltonian and any constraint. The only restriction is that \^d) has a positive overlap 
with the eigenstate (^f c \ of the constrained propagation. This restriction naturally implies a 
repulsive contribution to E g and an increase in the value of the energy. This algorithm can 
be made quite general and applied to a variety of interesting situations. 



In this sub-section we consider a simple numerical example illustrating the behavior of 
the various energy estimators. The particular example is the 2D Hubbard model on a 4 x 4 
lattice with 5 up spin and 5 down spin electrons. The exact ground-state energy of this small 




< a . 



(41) 



[(* T |* T )(0|0)]V2 
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system was obtained by direct diagonalization. For the intermediate coupling of U = 8t and 
t — 1, the energy is -17.51037.0 

We used the free-particle wave functions for both the constraint and the importance 
function in a series of CPMC calculations. As a variational wave function, the free-particle 
wave functions are quite inaccurate, yielding an energy of -11.50. We also used population 
sizes of 1000 to 2500 configurations, divided into two halves for independent left- and right- 
hand wave functions. Averages were computed over 30-100 blocks with a propagation time 
(number of steps times At) of 2 to 10 per block. We verified that we have reached the 
equilibrium state before computing averages and that the blocks were large enough to avoid 
difficulties with autocorrelations among individual energy estimates. All calculations were 
performed on single workstations, though extensions to large systems would require a parallel 
implementation. 

The various energy estimators are plotted as a function of the size of the time step At 
in Fig. 0. The exact ground-state energy is shown as a circle at the extrapolated At = 
limit. The two dashed curves illustrate the growth and mixed estimates E g and E m . E g is 
simply obtained from the change in overlap with iteration 



while E m is obtained by direct evaluation of Eq. [12]. Since the propagator is approximate, 
these two estimates coincide only in the limit of small At. 

As apparent from the figure, these two estimates lie slightly below the exact energy. The 
value of Eg, extrapolated to At = 0, is —17.517(2). During the course of this calculation, 
we also evaluated E&. In this case, Ed is small and positive, and adding this difference to 
Eg should produce a variational upper bound to the ground-state energy. Extrapolating to 
At = 0, we find E^ = 0.010(1), and hence E v = —17.506(2). The accuracy of the variational 
bound is quite surprising: the exact energy is recovered with an accuracy of two parts in 
10 4 , or better than 99.9% of the difference between the exact and trial state energies. 

We also plotted in Fig. |] the result of a direct calculation of E v . Since we have inde- 




(42) 
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pendent calculations of the left and right-hand states, it is possible to combine these into 
a direct calculation of the variational energy E v . In contrast to the other estimators, this 
calculation should yield a variational upper bound independent of the time step At. We 
find this to be true, but with a somewhat larger statistical error than the other estimators. 
There also appears to be some residual statistical bias resulting from the finite population 
size. This estimator may be more difficult to compute reliably for larger system sizes. 

The statistical errors for Ed also increase rapidly with system size. Primarily this is a 
result of a large statistical error in the denominator (Eq. |39| ). Particularly for our simple 
choices of trial states, the overlaps of the configurations representing the left- and right-hand 
population can vary dramatically. The alternative method of altering the constraint slightly 
(Eq. ^TJ) should produce a more favorable scaling with system size, though it remains to be 
demonstrated that this is practical for very large simulations. 

For the 4x4 case, we explicitly changed the constraint by introducing a finite value of 
a. We could achieve a variational upper bound by setting a = 0.0005; for a = 0.0003 
we could explicitly see that the sign of the overlap (^d^d) could lead to a violation of 
the upper bound. For a = 0.0005 and a time step of 0.005, we obtain a mixed estimate 
E m = —17.518(3) and a growth estimate E g = —17.505(3). Recall that only E g provides an 
upper bound in the limit of zero time step. However, the small difference E m — E g indicates 
the accuracy of the solution. Extrapolating to zero time step yields E g = —17.510(10). 

We also considered a 6x6 lattice with 13 spin up and 13 spin down electrons, again for 
U =8. We are unaware of any exact or QM C calculations for this system size at this filling. 
These larger system size results are meant to serve as guides for future use rather than 
exhaustive calculations. They were obtained on single-cpu workstation over the course of a 
few days. 

For the 6x6 system the constant a must be decreased significantly. This is rather natural 
as one would expect it to scale roughly with a small power of the number of single-particle 
orbitals. Again we use approximately 1000 configurations averaged over 30-100 individual 
blocks with a total propagation time of 2-4 per block. Here it is not clear if the original 
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choice of a = provides a variational upper bound, estimates of Ed bracket zero within the 
statistical errors of the calculation. For a = we obtain E m — E g — —36.05(05). 

Increasing the constant a to 10 -6 provides a variational upper bound. In this case we 
obtained E m = —35.75(05) and E g = —34.55(10) for a time step of At = 0.005. Extrap- 
olating to At = yields E m = -35.80(05) and E g = -35.25(20). It is possible that this 
bound could be further improved by using a somewhat smaller value of a. Again, the few 
per cent difference between E g and E m indicates the accuracy of the calculation. 

V. CONCLUDING REMARKS 

We presented several differences between the constrained-path and fixed-node Monte 
Carlo methods, some major and some minor. The most significant consequence of these 
differences is the mixed estimator in the CPMC method not being an upper bound to the 
exact energy as it is in the fixed-node method. Alternate ways of producing an upper bound 
have been introduced. 

While not an upper bound, the mixed estimator in the CPMC method was argued 
to be very near the exact answer. Experience shows it is almost always above the exact 
answer, and in cases where the CPMC results have been compared to fixed-node results, the 
CPMC mixed estimates of the energy always lie closer to the exact answer than the upper 
bound produced by the fixed-node method. El Presumably this accuracy is a consequence 
of the quality of the estimates of the wave function. As a rule of thumb, we find that the 
fewer nodal crossings the more accurate is the prediction of the energy. This observation 
is supported by the method discussed in Section IV to correct the mixed estimate so it 
produces an upper bound. This method depends on computing a contribution from those 
walkers thrown away; if none are thrown away (after the sampling is from the limiting 
distribution), then the CPMC method in fact becomes exact. 

There are several other differences between the methods worth mentioning. In one con- 
tinuous spatial dimension, coincident planes (x, = Xj for all i and j corresponding to the 
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same spin species) exhaust the nodal surface set for local potentials V; therefore, in the 
fixed-node method one can get the exact ground-state energy by using any nodal surface 
with this property. The free-fermion wave function suffices. For general lattice fermion 
problems, even in one spatial dimension, the situation is more complicated: there are extra 
nodal surfaces which are not coincident planes. For certain classes of Hubbard-like models 
the latter exhaust the whole nodal set; therefore it is possible to avoid the sign problem and 
get the exact solution.0 On the other hand, the lattice version of the fixed-node method0 
always provides a variational upper bound to the exact ground-state energy regardless of 
the dimensionality, unless the constraining state is the exact ground-state wave function. 
We have been unable to develop a similar understanding for the constrained-path method 
where in one dimension we observe an absence of nodal crossings and the mixed estimate 
of the energy agreeing with exact results to statistical accuracy.0 We comment that care 
must be taken in using nodal crossing rule of thumb. The accuracy of the Trotter approxi- 
mation is controlled by the size of At. If At is large, the approximation is poor, and nodal 
crossings can be induced into a problem for which there is no sign problem. If At is too 
small, the propagation through phase space is too slow. For a poor choice of importance 
sampling population control can sweep away walkers that should cross the surface before 
they actually do. 

We also remark that in the CPMC method there is much less need to perform a variational 
optimization of \^t) through Jastrow or Gutzwiller factors as seems to be necessary in the 
fixed-node method.! This optimization process does not affect the nodal surface but does 
reduce the energy of the starting configuration. While there is always some advantage in 
doing this, the results of the CPMC method display considerable robustness to the choice 
of the constraining wave function which is also typically used as the starting configuration. 
Simple choices, like free-fermion wave functions, seem to work well. Quite different choices 
of I^t) usually give satisfyingly similar results.!! 

In closing, we remark that some of our observations about the mixed estimator for the 
energy might be useful in constructing and assessing estimation procedures used in the 

20 



LA-UR-98-5289 



February 1, 2008 



submitted to Phys. Rev. B 



standard auxiliary-field projector quantum Monte Carlo (AFQMC) method. In that 
method, the energy or more generally some observable O is typically estimatedlMl'Ill from 




and then for large r this expression is either averaged over several values of r' or evaluated 
at just t' = t/2. Clearly, the latter procedure may be preferable, even though the former 
may have lower variance, as this estimator can be rewritten as 



revealing that for r' 7^ r it is basically just a mixed estimator. For estimating the energy 
or an observable that commutes with H, the utility of this estimator depends on how close 
either or \^r) approaches the ground state wave function. For estimation of observables 
that do not commute with the Hamiltonian H, its utility depends on how close both 
and \^r) approach the ground state wave function. 



We gratefully acknowledge illuminating discussions with Erik Koch on the variational 
aspects of the constrained-path method. Most of this work was supported by the Department 
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and the Research Corporation. 



We now introduce an exactly solvable fermion model to illustrate and visualize several 
features of the constrained-path method. This model has the Hamiltonian 



0(t, t') 



(*l\Q\Vr) 
(*l|*h> 



(44) 
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APPENDIX: ILLUSTRATIVE EXAMPLE 



H = -t (c\c 2 + C2C3 + e 2 c x + c\c 2 ) + U (n x n 2 + ^2^3) = T + V , 



(Al) 
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and corresponds to spinless fermions coupled through a nearest neighbor repulsive interaction 
[Zona three site lattice with open boundary conditions. In the following we will concentrate 
on the two-particle solutions for which the ground state is 

1 



l*o> 



2t 2 + {U- E y 



t {c{ c \ + 44) 



(U-E )c\ct |0) 



(A2) 



and E , the ground-state energy, equals U/2 — J(U /2) 2 + 2t 2 . 

Since we want to study explicitly the (imaginary-time) evolution of the distribution of 
Slater determinants which arise as a consequence of the constrained-path approach, we need 
to have some way of parameterizing the differentiable manifold of Slater determinants of 
two particles that has dimension N(M — N) = 2(3 — 2) = 2. We can parameterize a state in 
the two-particle Hilbert space Ti.2 that belongs to the set of normalized Slater determinants 
by the 3x2 matrix 



cos(#i - 6 2 ) 

COS 6\ COS 6*2 

sin 6*i sin 6 2 



(A3) 



Therefore, the two angles (9\, 82) specify a point in the manifold, and any state belonging 
to the Hilbert space H.2 having support in that manifold can be represented by \l/[0i, 62] = 
(01^). In general, the usual property of a fermion wave function to be totally antisymmetric 
in spin-coordinate space is lost in this new representation. For instance, the ground state of 
our model fermion system is given by 

1 



#o[0i,0 2 ] 



2t 2 



(U - E ) 



[cos(0i - 2 ) (t cos 2 + (17 - Eq) sin 6 2 ) + t sin(0 2 - 6 X )\ 



(A4) 



which is neither symmetric nor anti-symmetric under permutations of 9\ and 9 2 - Contour 
plots of it can be found in Fig. |3| for different values of the interaction strength U. In the 
non-interacting case, U = 0, the ground state is a single Slater determinant represented by 
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the matrix 



corresponding to the point 0\p = @2,t — 7r/6, #2,t = arccos 1/V3 in the manifold. 

For a given value of U we would like to project out the ground state from this non- 
interacting state. For simplicity, we will leave out importance sampling and use a first order 
Trotter decomposition for the short time propagator 



-AtH 



10) = e 



-AtV „-AtT\ 



+ 0{Ar 2 ) , 



(A6) 



and stochastically iterate this expression for each walker. In a matrix representation this it- 
eration is equivalent to a product of non-commuting random matrices, i.e. $(r) = U(t)<& t = 
n™=i (exp[— AtV] exp[— ArT]) $t, where the symmetric matrix exp[— ArT] is given by 

u + 1 v u — 1 1 



-AtT 



V2 



2 

72 " V2 
u — 1 t> w + 1 



(A7) 



V 2 v/2 2 7 

with u = cosh(v / 2 tAr), t> 2 = u 2 — 1. After the use of the discrete 
([/ > 0) Hubbard- Stratonovich transformation exp[— ArUrijUj+i] = \ exp[— ArU{rij + 
rij+i)/2] J2 x =±i exp[/ix(rij — J^+i)], the interaction part of the propagator at any imaginary- 
time slice % is represented by one of the 4 diagonal random matrices (each chosen with 
probability 1/4, i.e. P(x) = 1/4) 
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where cosh/x = exp[Ar{7/2]. Although U{j) is a random matrix, its determinant det U(t) = 
exp[— 2tII] is not a random number. 

After a short propagation (Eq. |A6| ) , each point of the Slater determinant manifold, 
representing the state of the system, performs a Brownian walk in 9— space. Therefore, one 
can consider 9\{r) and 9 2 (t) (such that #i(0) = 9\ t T, and #2(0) = 02,t) as a se t of random 
variables defining a random walker in imaginary time. However, in the CPMC method, to 
avoid the fermion sign problem, we constrain the walker with a constraining state which 
here we choose it to be \^fo(U = 0)). In other words, each time Eq. [AH is iterated, we ask, 
Has the sign of det $£<I>(t) changed ? If it has, then we kill the walker. 

For the present example, the determinant can be computed analytically: At any time 
step n > 1 



det 



"i 11 ~\~ 1) 

®t®(t)} = Mi f 2 + 92 + 2«i7i M , (A9) 



where the functions gi and hi are most simply obtained from the backwards recursion 
relations 

fi = aiPi (u + 1) f i+ i + frji (u - 1) g i+1 + lacfi v h i+1 

gi = ai(3i (u - 1) f i+1 + /3iji (u + 1) g i+ i + 2a^i v h i+1 (A10) 

hi = aiPi v /j+i + v g i+1 + lacii u K+i 

with f n+ i = g n+ i = h n+ i = 1 being the initial conditions. 

Clearly the determinant Eq. is always positive, independent of the values of n, At and 
U, which means that there is no sign problem. Nevertheless, this example helps to illustrate 
several important issues. The first is that the exact ground-state energy can be stochastically 
obtained (as we will show below) even if the nodal surface of is approximate. The 

second issue deals with the practical numerical implementation of the method. A "false" sign 
problem (i.e. det $£$(r) < 0) can occur as a consequence of numerical round-off errors. 
In fact, the use of matrix stabilization techniques is crucial to avoid such phenomenon.0 
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Figure f| shows the energy mixed estimator E m as a function of imaginary time. All 
walkers start (r = 0) at the point (0i j t,&2,t) m #-space. That means that E m (r = 0) = 
U/2 — V2t, and as time evolves, E m — > Eq if the constraint is not evoked. As mentioned 
above this is the case if matrix stabilization techniques are properly used, otherwise we 
observe "false" nodal crossings. A typical random walker 6(t) is displayed in Fig. |5]. This 
figure clearly denotes that most of the time the walker prefers to stay near the upper right 
corner of the space. This behavior is evidenced in the lower panel of the same figure, where 
it is shown that when the walker is at the upper right corner in #-space the overlap with the 
exact ground state is maximum (in fact, it is almost one). 
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FIGURES 

FIG. 1. Contour plots of the two-fermion ground-state wave function (^l^o) for various values 
of the interaction strength U/t. Notice the differences among the nodal structures. 

FIG. 2. Energy estimators as a function of time step for the 4x4 Hubbard model described 
in the text. Solid symbols indicate estimators which are variational upper bounds to the exact 
energy. 

FIG. 3. Contour plots of the two-fermion ground-state wave function (^l^o) f° r various values 
of the interaction strength U/t. 

FIG. 4. Energy mixed estimator as a function of imaginary time averaged over N w walkers. 
The horizontal dashed line indicates the value of the exact ground-state energy Eq = — 0.19615t. 

FIG. 5. The upper panel shows a typical random walk in 6>-space. At r=0 the walker starts at 
(0i,Tj 02,t)- The random walk never crosses the nodal surface of (4>\^t)- The lower panel displays 
the overlap of the walker with the exact ground state, and the distance of the walker from the 
origin. Notice the clear correlation between these two quantities. 
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